rm(list=ls())

require(RColorBrewer)
require(logistf)
require(data.table)
require(MASS)
require(glm.predict)
require(kitagawa)
require(magicaxis)
require(tidyverse)
require(lattice)
require(latticeExtra)
require(grid)

##############################################
#Source functions necessary for calculations #
##############################################
source("functions.R")

##############################################################
#Read UN population data for S. Korea                        #
#Reported in 1000s                                           #
#Source: https://population.un.org/wpp/Download/Standard/CSV/#
##############################################################
world.population <- read.csv("~/Dropbox/mortality/covid19/data/WPP2019_PopulationByAgeSex_Medium.csv")

###################################################################################################################################################
#S. Korean ASDR by cause of death
#Source: #
###################################################################################################################################################
cod <- read.csv("~/Dropbox/mortality/covid19/data/cod_kor.csv")
cod <- c(as.numeric(subset(cod, Year==2016)[-1]))
names(cod) <- c("Cancer", "Heart Disease", "Cerebrovascular Disease", "Pneumonia", "Intentional Self Harm", "Diabetes Mellitus", "Chronic Lower Resp. Tract Disease", "Liver Disease", "Traffic Accident", "Hypertensive Disease")

#########################################################################
#South Korea COVID-19 deaths by age group                               #
#Source: https://www.cdc.go.kr/board/board.es?mid=a30402000000&bid=0030 #
#########################################################################
starting.ages <- seq(0,80,10)
population <- population.fxn(world.population, "Republic of Korea", 2020, starting.ages)
deaths <- c(0,0,0,2,3,15,34,70,113)
cases <- c(138,582,2928,1139,1413,1951,1343,706,483)
population.standard <- prop.table(population.fxn(world.population, "Republic of Korea", 2016, starting.ages))
start.date <- "2020-01-21"
end.date <- "2020-04-21"

mortality.outcomes <- covid19.mortality.rate.fxn(starting.ages, population, cases, deaths, case.fatality.rates, start.date, end.date)
mx <- mortality.outcomes$mx.firth
graph.mx.fxn(mx)
create.table.mx.fxn(mx)
graph.alternative.mx.fxn(mortality.outcomes$mx.firth, mortality.outcomes$mx.nb, mortality.outcomes$mx.poisson)
elasticity <- mortality.outcomes$elasticity
create.table.elasticity.fxn(elasticity)
deaths.data <- mortality.outcomes$deaths.data.wide

prevalence.max <- 5*200210/sum(population.fxn(world.population, "Spain", "2020", starting.ages)) #Number of confirmed cases in Spain was 200,210 on 21 April 2020 (WHO Sit Rep 92)
mortality.undercount.scenario <- seq(0,1,length.out=100)
asdr <- asdr.mat.fxn(prevalence.max, mortality.undercount.scenario, mortality.outcomes, population)
graph.asdr(asdr)